Background

The body condition and growth of Eastern Baltic cod (Gadus morhua) has declined steadily since the regime shift in the early 1990’s to a degree that the stock now can be viewed as collapsed. Several hypotheses have been put forward, including changes in overlap with pelagic prey (e.g. Gårdmark et al, 2015), reduced oxygen levels decreasing habitat quality and leading to contraction of the distributional range thus increasing competition (e.g. Casini et al, 2016), increased competition for benthic food sources with flounder (Orio et al, 2019; Orio 2020) as well as increased intraspecific competition and growth bottlenecks within the population (Svedäng & Hornborg, 2014).

However, these potential explanatory variables have not been evaluated on a fine spatial scale, even though factors such as competition, food availability and local habitat quality likely occur on fine spatial scales. Instead, averages over larger spatial areas (e.g. ICES subdivisions) have been the variables in previous comparisons. Moreover, the ability of each of proposed explanatory variable (linked to a hypothesis) to explain variation in the condition of cod has not been compared in a standardized way, and has not been contrasted to residual spatial and spatiotemporal variation.

Aim

In this study I have compiled data for individual-level condition and matched that with predictor variables (on a haul level) representing different ecological hypotheses regarding drivers of variation in cod condition.

To account for spatial and temporal autocorrelation using data at this scale, I in apply spatiotemporal predictive-process GLMMs using the R-package sdmTMB. This modeling framework allows evaluation of how much of the variation in condition can be explained by covariates, spatial (unmeasured variation in condition that is stable over time) and spatiotemporal variation (unmeasured variation in condition that changes between years).

In this script, I aim to find a good default model explaining spatiotemporal variation in cod condition, without any of the suggested explanatory variables (the rest will be done in a different script). Then I show examples of two models with additional explanatory variable. The project currently lives here. Details follow below.

Methods

Modeling framwork

In fishes, weight is typically assumed to vary log-normally around an average allometric function of length: \(w=al^b\), where \(w\) is weight in grams, \(l\) is length in cm, \(b\) is the allometric length exponent and \(a\) is the condition factor in unit \(g/L^b\) (Froese et al., 2014; Grüss et al., 2020). Typically this relationship is linearized by taking logs on both sides: \(\operatorname{log}(w)=a+b\operatorname{log}(l)\). Le Cren’s condition index is defined as the residuals from this length-weight relationship.

We model this individual-level relationship with a spatiotemporal GLMM of the form:

\[ \operatorname{log}(w_{s,t}) = \alpha_0 + \alpha_mS_m + \alpha_uS_u + \omega_s + \epsilon_{s,t} + \sum^{n_k}_{k=1}\gamma_k\boldsymbol{X_k} + \beta \operatorname{log}(l) \]

where \(\alpha_0\) is the female condition factor, \(\alpha_m\) is the offset for males and \(\alpha_u\) for unsexed individuals. \(\omega_s\) and \(\epsilon_{s,t}\) represent spatial and spatiotemporal random effects, respectively. \(\boldsymbol{X_k}\) is a matrix of \(n_k\) measured additional covariates and \(\gamma_k\) is the effect of the \(k\)-th additional covariate. \(\beta\) is the length-coefficient, corresponding to the allometric exponent. The spatial and spatiotemporal random effects are assumed to be drawn from a multivariate normal distribution:

\[ \omega \sim \operatorname{MVNormal}(0, \Sigma_\omega)\\ \epsilon \sim \operatorname{MVNormal}(0, \Sigma_\epsilon) \]

I also consider the spatiotemporal random effects to be drawn from a multivariate normal distribution following an AR1 process:

\[ \delta_{t=1} \sim \operatorname{MVNormal}(0, \Sigma_\epsilon)\\ \delta_{t>1} = \phi\delta_{t-1} + \sqrt{1-\phi^2}\epsilon_t, \epsilon_t \sim \operatorname{MVNormal}(0, \Sigma_\epsilon) \]

In the spatial and spatiotemporal random fields, \(\Sigma_\omega\) and \(\Sigma_\epsilon\) are covariance matricies, where the covariance (\(\Phi(s, s')\)) between spatial points \(s\) and \(s'\) is given by a Matérn function:

\[ \Phi(s, s') = \tau^2/\Gamma(\nu)2^{\nu-1}(\kappa d_{jk})^{\nu}K_\nu(\kappa d_{jk}) \] where \(\tau^2\) is the spatial (marginal) variance.

This model (first equation) can be viewed as an approximation of Le Cren’s condition index (Grüss et al., 2020), as the log of the condition factor, i.e. \(\operatorname{log}(a)\) or the constant in the allometric relationship, can be defined as: \(\operatorname{log}(a) = \alpha_0 + \alpha_sS + \omega_s + \delta_{s,t} + \sum^{n_k}_{k=1}\gamma_k\boldsymbol{X_k}\). Thus, Eq. 1 is a model for a spatially and temporally varying condition factor.

Finding a default model

I aim to find an appropriate model structure for Eq. 1, to which I later can include and evaluate additional covariates (\(\gamma_k\boldsymbol{X_k}\)). To find this default model, I fit a model with only length at sex as covariates, and assume errors are Student-t distributed to account for the presence of heavy tails.

I fit different versions of this model:

  • Fixed effect of year
  • Time-varying intercept through random walk
  • Spatial trend, which is a random field representing a slope over the time variable.

These models are fit with the spatiotemporal field being independent each year or following an AR1 process.

Data

Individual length and weight

The data are individual-level measurements of length and weight of Baltic cod in quarter 4 between the year 1991-2019, and stem from the Baltic International Trawl Survey (BITS) (can be downloaded from DATRAS). These data are cleaned and merged with additional covariate data in this script.

Haul-level covariates

The covariates I currently consider are:

  • Density of cod (by length group, currently >30cm, ≤30 cm or total). Abundance CPUE available for each haul from the BITS survey.
  • Density of flounder (by length group, currently >20cm, ≤20 cm or total). Abundance CPUE available for each haul from the BITS survey.
  • Oxygen concentration. Possible to link ICES oceanographic data data to nearest haul in BITS
    • added, but very coarse resolution. Have asked the Swedish Meteorological and Hydrological Institute for model output since current samples are sometimes too far from the haul location
  • Sprat CPUE. These data stem from the Baltic International Acoustic Survey (BIAS).
  • Herring CPUE. These data stem from the Baltic International Acoustic Survey (BIAS).
    • Each haul in the condition data is assigned with an abundance of sprat and herring based on the ICES rectangle the haul was in, because we don’t have sprat and herring abundance on a finer scale than rectangle. Note also these data are not really available online, because some calculations have been made to get abundances per rectangle…
  • Maybe: temperature as a base-covariate through its effect on growth, metabolism, digestion etc.

The covariates are standardized to have a mean of 0 and a standard deviation of 1, to allow for direct comparison with the spatial and spatiotemporal standard deviation (Thorson, 2015; Grüss et al., 2020)

Explore data

library(tidyverse); theme_set(theme_classic())
library(tidylog)
library(viridis)
library(sdmTMB)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)

# For adding maps to plots
world <- ne_countries(scale = "medium", returnclass = "sf")

Now read data:

d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/clean_for_analysis/mdat_cond.csv")

# These should be factors...
d <- d %>% mutate(year_f = as.factor(year), # Actually I'm modelling this as a continous variable for now
                  sex_f = as.factor(sex))

# Since I also use year as a continuous variable, I will center it to start at 1991 (start of time series)
# (*** Not used currently, I opted for the factor effect of year instead as that lead to better convergence as in smaller max gradient. To use a continuous year variable I had to use fewer knots, and then the models without year had some issues with convergence (m3 below))
d <- d %>% mutate(year_ct = year - 1991)

Read the prediction grid:

pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/clean_for_analysis/pred_grid.csv")

# SA: You likely want to work in UTMs so distance is constant with space
# ML: I did not think so much about this, but that sounds like a good idea. Although for now I can't choose which of the two UTM zones to use... (I did change from X and Y to lon and lat to make it more clear in case I switch
# later!)
pred_grid <- pred_grid %>%
  mutate(ln_length_cm = log(1),
         year_f = as.factor(year), # Using year as a continious variable in the model as there seems to be a
         # linear decline
         sex_f = "F") %>% # For now we'll predict changes in the intercept ("female condition factor")
  mutate(X = lon,
         Y = lat,
         year_ct = year - 1991) %>% 
  filter(year %in% c(unique(d$year)))

We can now plot the Fulton K condition factor to get a glimpse of what we might expect (but in the main model we use Le Cren’s condition index).

# Plot "Fulton K" in space and time 
# It's a little bit difficult to see because we have many extreme observations
ggplot(d, aes(Fulton_K)) + geom_histogram()
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


# Plot "Fulton K" in space and time 
# It's a little bit difficult to see because we have many extreme (high condition index) observations
ggplot(d, aes(lon, lat, color = Fulton_K)) + 
  geom_point(size = 1.5, alpha = 0.8) + 
  facet_wrap(~ year, ncol = 6) +
  scale_color_gradient2(midpoint = mean(d$Fulton_K)) +
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58))

Ok, there’s a clear temporal development in condition and the spatial coverage of data varies by year (fewer data in the beginning of the time series).

Find the default model

Compare models

After some exploration of all models considered here, I found that ~110 knots was the highest I could and still get convergence. I use run_extra_optimization (with various settings) for all models to reduce max gradient, but with higher # of knots I get an error I do not know how to resolve except with fewer knots:

Warning messages: 1: In sqrt(diag(cov)) : NaNs produced 2: The model may not have converged: non-positive-definite Hessian matrix.

# SA: good for now; could increase knots eventually:
# ML: Thanks! I did increase from 70 up to 110, but see above 

spde <- make_mesh(data = d, xy_cols = c("lon", "lat"), n_knots = 110,
                  type = c("kmeans", "cutoff", "cutoff_search"))
plot(spde)

I assume use a Student-t distribution with 3 degrees of freedom (default), because initial exploration suggest presence of heavy tails in residuals. Now, fit the models outlined above.

# SA: note that the degrees of freedom on the Student-t are fixed out 3 currently; we can tweak that so you could specify it if you want.
# ML: Yes, I saw that! To me at least it seems that the Student model drastically improves fit by accounting the the heavy tails, but maybe assumes too heavy tails? I fit a simple glm (weight ~ length) in brms and estimate nu to be around 5.7.
# SA: working in log space and Student-t might be fine here, though it's possible that something like Gamma(link = "log") or lognormal might capture the mean-variance relationship better and not require the transformation.
# ML: Ok, I see! In this case the main reason for taking the logs is transform the allometric model to a linear form, so maybe that wouldn't work

# Model 1
m1 <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f + year_f, data = d, time = "year",
             spde = spde, family = student(link = "identity"), ar1_fields = TRUE,
             include_spatial = TRUE, spatial_trend = FALSE, spatial_only = FALSE,
             silent = TRUE) 
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.405740254218188.

m1b <- run_extra_optimization(m1, nlminb_loops = 1, newton_steps = 2)
max(m1$gradients)
#> [1] 0.4057403
max(m1b$gradients)
#> [1] 4.154117e-10

# Model 1.1
m1.1 <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f + year_f, data = d, time = "year",
               spde = spde, family = student(link = "identity"), ar1_fields = FALSE,
               include_spatial = TRUE, spatial_trend = FALSE, spatial_only = FALSE,
               silent = TRUE) 
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.57370793938685.

m1.1b <- run_extra_optimization(m1.1, nlminb_loops = 0, newton_steps = 1)
max(m1.1$gradients)
#> [1] 0.5737079
max(m1.1b$gradients)
#> [1] 7.959633e-06

# Model 2
# m2 <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f, time_varying = ~ 1 + sex_f,
#              data = d, time = "year",
#              spde = spde, family = student(link = "identity"), ar1_fields = TRUE,
#              include_spatial = TRUE, spatial_trend = FALSE, spatial_only = FALSE,
#              silent = TRUE)

# Model 3
m3 <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f, data = d, time = "year",
             spde = spde, family = student(link = "identity"), ar1_fields = TRUE,
             include_spatial = TRUE, spatial_trend = TRUE, spatial_only = FALSE,
             silent = TRUE) 
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.518019836650495.

m3b <- run_extra_optimization(m3, nlminb_loops = 0, newton_steps = 1)
max(m3$gradients)
#> [1] 0.5180198
max(m3b$gradients)
#> [1] 1.352712e-06

# Model 3.1 - spatial trend instead of AR1 spatiotemporal
m3.1 <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f, data = d, time = "year",
               spde = spde, family = student(link = "identity"), ar1_fields = FALSE,
               include_spatial = TRUE, spatial_trend = TRUE, spatial_only = FALSE,
               silent = TRUE) 
#> Warning: The model may not have converged. Maximum final gradient:
#> 1.17092374178321.

m3.1b <- run_extra_optimization(m3.1, nlminb_loops = 0, newton_steps = 1)
max(m3.1$gradients)
#> [1] 1.170924
max(m3.1b$gradients)
#> [1] 2.708515e-05

Regarding model 2 here: I’m a bit uncertain on what model I’m actually fitting here. Since I have the sex effect, which represent the intercept when sex="F" and an offset when "M" or "U", I’m guessing it would only make sense to have all those coefficients as time varying. It doesn’t converge, however, and I’ve also tried: ar1_fields = FALSE. I get:

Warning message: The model may not have converged: non-positive-definite Hessian matrix.

I then tried omitting the sex-effect and modeled only the intercept as time-varying:formula = ln_weight_g ~ ln_length_cm, time_varying = ~ 1. But here I get: Error in solve.default(h, g) : system is computationally singular: reciprocal condition number = 2.12999e-16

I’m proceeding with m1b, m1.1b, m3b and m3.1b though. Look at the residuals:

df <- d

df$residuals_m1b <- residuals(m1b)
df$residuals_m1.1b <- residuals(m1.1b)
df$residuals_m3b <- residuals(m3b)
df$residuals_m3.1b <- residuals(m3.1b)

qqnorm(df$residuals_m1b); abline(a = 0, b = 1)

qqnorm(df$residuals_m1.1b); abline(a = 0, b = 1)

qqnorm(df$residuals_m3b); abline(a = 0, b = 1)

qqnorm(df$residuals_m3.1b); abline(a = 0, b = 1)

While the Student-t models looks a lot better (than the Gaussian [not shown]) they could perhaps be improved further by tweaking the df parameter (see comment above). Otherwise, they look identical to me (I am a bit surprised by that actually). I will next compare their AIC:

AIC(m1b)
#> [1] -169633.5
AIC(m1.1b)
#> [1] -169626.3
AIC(m3b)
#> [1] -169544.8
AIC(m3.1b)
#> [1] -169546.2

Model m1b, i.e. with an autoregressive spatiotemporal field and with a fixed (factor) year effect is most parsimonious and seems like a good option to proceed with now (given that I can’t detect any different in the QQ plot between the alternative models).

Parameter estimates

We can also look at the AR1 parameter, and sex-coefficients, to ensure they are warranted.

# SA: Note that this is in -Inf to Inf space, although you are probably aware of that; `2 * plogis(ar_phi) - 1` will transform it
# SA: also note that you can use (1) the already calculated sdreport and (2) the as.list S3 method for a handy named list... no more row.names()!:
# ML: Right, that's what you wrote in the e-mail, forgot to correct that.
# ML: Perfect, will use as.list on the report. Thanks!
m1b_est <- as.list(m1b$sd_report, "Estimate")
m1b_se <- as.list(m1b$sd_report, "Std. Error")
# Transform back to -1 to 1 scale
2 * plogis(m1b_est$ar1_phi) - 1
#> [1] 0.162155
2 * plogis(m1b_est$ar1_phi + c(-2, 2) * m1b_se$ar1_phi) - 1
#> [1] 0.05511056 0.26551430

Strong support for it (and the AIC-selected model also had it). What about the sex_f-effect?

# Females
m1b_est$b_j[1]
#> [1] -4.431723
m1b_est$b_j[1] + c(-2, 2) * m1b_se$b_j[1]
#> [1] -4.470238 -4.393209

# Males
m1b_est$b_j[3]
#> [1] -0.005740679
m1b_est$b_j[3] + c(-2, 2) * m1b_se$b_j[3]
#> [1] -0.006983113 -0.004498246

# Unsexed
m1b_est$b_j[4]
#> [1] -0.08561291
m1b_est$b_j[4] + c(-2, 2) * m1b_se$b_j[4]
#> [1] -0.09032741 -0.08089842

The condition appears different between sexes. The average condition for a female is -4.43, and for males it is -4.44. Given the length-coefficient (\(\beta\) or m1b_est$b_j[2]), we can exemplify the difference in weight:

# How much does the weight of a 30 cm female and male differ?
exp(m1b_est$b_j[1])*30^m1b_est$b_j[2]
#> [1] 282.6674
exp(m1b_est$b_j[1] + m1b_est$b_j[3])*30^m1b_est$b_j[2]
#> [1] 281.0494
# Females are approximately 0.5 g heavier, not that much maybe [in other models ~1 g for the same size]! I could consider perhaps dropping that term and then compare that model using AIC.

Plot the residuals

First we can plot them against predictions and then over length:

pred_m1b <- predict(m1b)
df$pred_m1b <- pred_m1b$est

ggplot(df, aes(pred_m1b, residuals_m1b)) +
  geom_point(shape = 21, fill = "black", color = "white")


ggplot(df, aes(ln_length_cm, residuals_m1b)) +
  geom_point(shape = 21, fill = "black", color = "white")

They look a little bit funny at small lengths (and small predicted weights). I think that’s because a) the small ones are quite rare and b) data are in cm size-classes (a bit coarse for a ~5 cm fish, hence many are assigned the same length). Besides that I do not think they are too troublesome.

Check the residuals for the Student-t model on a map.

ggplot(df, aes(lon, lat, colour = residuals_m1b)) +
  geom_point(size = 1) +
  facet_wrap(~year, ncol = 6) +
  scale_color_gradient2() +
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58))

Maybe some clustering remains… But overall OK I think!

Plot predictions on map

Now we can predict and plot estimates using all fixed and random effects on pre-made grid. This grid is created by doing an expand.grid over the survey range, then filtering out areas that are actually in the ocean using ICES shapefiles. Lastly some areas are too deep for sampling (-135 m). I’ve added a depth column so that I can make those predictions NA so it’s clear they are different from e.g. land and islands (and so that the color gradient isn’t going too far because of the extreme predictions at these depths).

p <- predict(m1b, newdata = pred_grid)

# Replace too-deep predictions with NA
p <- p %>% mutate(est2 = ifelse(depth < -130, NA, est), # prediction (fixed + random)
                  eps_st2 = ifelse(depth < -130, NA, epsilon_st), # spatiotemporal effects
                  omega_s2 = ifelse(depth < -130, NA, omega_s), # spatial random effect
                  est_non_rf2 = ifelse(depth < -130, NA, est_non_rf), # fixed effects
                  zeta_s2 = ifelse(depth < -130, NA, zeta_s)) # spatial trend

Plot the predicted condition with fixed and random effects:

ggplot(p, aes(lon, lat, fill = est2)) +
  geom_raster() +
  facet_wrap(~year, ncol = 6) +
  scale_fill_viridis(option = "magma", 
                     name = "log(condition factor)") + 
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58)) +
  ggtitle("Prediction (random + fixed)")

Plot the spatiotemporal random effects:

SA: Note that these are the spatial (omega) + spatiotemporal (epsilon) effects. If you just want the spatiotemporal ones, then you want p$epsilon_st. I have done that above. ML: Oh didn’t realize that. Good, perhaps that was the reason they seemed similar… SA: You likely want a diverging colour scheme for these since they will be centered on 0. I have done that. ML: Yes makes sense and looks good! SA: typically each time slice would be centered on zero, however you have a moderately large AR1 phi parameter, which is letting it wander off a bit SA: an alternative might be the spatial_trend (since what you have look roughly linear) likely instead of the AR1 form, I would be interested to hear if that works well. If everything is generally decreasing, you could also try just having a fixed effect for year, which might be simplest of all and useful for inference… or let the intercept be a random walk through time: formula = ln_weight_g ~ 0 + ln_length_cm, time_varying = ~ 1 SA: hopefully that’s right, I haven’t tried it here SA: All depends on what the end goal is for inference and what is most useful. ML: Those are good ideas! I have implemented that above.

ggplot(p, aes(lon, lat, fill = eps_st2)) +
  geom_raster() +
  facet_wrap(~year, ncol = 6) +
  scale_fill_gradient2(name = "eps_st") + 
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58)) +
  ggtitle("Spatiotemporal random effects")

Plot the spatial random effects:

ggplot(filter(p, year == 2000), aes(lon, lat, fill = omega_s2)) +
  geom_raster() +
  scale_fill_gradient2(name = "omega_s") + 
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58)) +
  geom_point(data = d, aes(lon, lat), color = "black", inherit.aes = FALSE, size = 0.1) +
  ggtitle("Spatial random effects + data")
#> filter: removed 119,140 rows (97%), 4,255 rows remaining

The spatial random field seems to largely follow depth, which makes sense since it’s constant across years and reflects e.g. oxygen concentration and overall habitat quality (deep areas are even anoxic).

baltic_sea <- getNOAA.bathy(lon1 = min(d$lon), lon2 = max(d$lon),
                            lat1 = min(d$lat), lat2 = max(d$lat), resolution = 15)
#> Querying NOAA database ...
#> This may take seconds to minutes, depending on grid size
#> Building bathy matrix ...

autoplot(baltic_sea, geom = c("r", "c")) +
  scale_fill_gradient2(low = "darkblue", high = "gray", midpoint = 0)

SA: there are also the p$est_non_rf for just the fixed effects (and time-varying effects) if you want to look at that ML: Sounds good, will plot that too!

Plot fixed effects on grid (note this is for female cod). If using year as a continuous variable, the overall result is not that different.

# Filter 1 year if I end up with a model like this where all years are identical
ggplot(p, aes(lon, lat, fill = est_non_rf2)) +
  geom_raster() +
  facet_wrap(~year, ncol = 6) +
  scale_fill_viridis(name = "log(condition factor)") + 
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58)) +
  ggtitle("Fixed effects")

# A bit clumpsy way to append the intercept (year 1991, female) to the other year-effects (coefficients 5:32)
year_intercept <- data.frame(Year = sort(unique(d$year)),
                             "est" = c(as.list(m1b$sd_report, "Estimate")$b_j[1],
                                       as.list(m1b$sd_report, "Estimate")$b_j[5:32] + as.list(m1b$sd_report, "Estimate")$b_j[1]),
                             "std.error" = c(as.list(m1b$sd_report, "Std. Error")$b_j[1], 
                                             as.list(m1b$sd_report, "Std. Error")$b_j[5:32] + as.list(m1b$sd_report, "Std. Error")$b_j[1]))

year_intercept %>% 
  mutate(ymin = est - 1.96*std.error,
         ymax = est + 1.96*std.error) %>% 
  ggplot(., aes(Year, est)) + 
  ylab("Average log(female condition)") + 
  geom_point() + 
  geom_errorbar(aes(ymin = ymin, ymax = ymax), width = 0.4)
#> mutate: new variable 'ymin' (double) with 29 unique values and 0% NA
#>         new variable 'ymax' (double) with 29 unique values and 0% NA

Here we see in more detail how the average condition for a female has changed. 1996 was a particularly good year for condition, and it happens to also be with with relatively good oxygen conditions. After I saw this plot, I also tried fitting a model with a spline for year (centered to 0 at 1991), but it did not have a lower AIC than this model (but lower then m3b and m3.1b)

Now we want to refit the same model with the additional fixed effects outlined above.

Add covariates to default model

Here is an example of how the importance of additional covariates can be evaluated. To keep it simple I contrast only two covariates (cod and oxygen). The rest of this analysis (i.e. with all covariates) will be done in different rmarkdown file.

# Fit model with cod cpue as covariate
mcod <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f + cpue_cod_st + year_f, data = d,
               time = "year",
               spde = spde, family = student(link = "identity"), 
               ar1_fields = TRUE,
               include_spatial = TRUE, 
               spatial_trend = FALSE, 
               spatial_only = FALSE,
               silent = TRUE)
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.204543740353301.

# Run extra optimization steps to help convergence:
mcodb <- run_extra_optimization(mcod, nlminb_loops = 0, newton_steps = 1)

max(mcod$gradients)
#> [1] 0.2045437
max(mcodb$gradients)
#> [1] 5.647645e-06

# SA: interesting; you might try more knots (or sometimes fewer) or even a different seed on the make_spde function.
# ML I landed on fewer knots because I had convergence issues also with the first models with more knots (that I couldn't resolve with ´run_extra_optimization´, as with this model originally!)

# And one with oxygen concentration as covariate
moxy <- sdmTMB(formula = ln_weight_g ~ ln_length_cm + sex_f + mean_DOXY_st + year_f, data = d,
               time = "year",
               spde = spde, family = student(link = "identity"), 
               ar1_fields = TRUE,
               include_spatial = TRUE, 
               spatial_trend = FALSE, 
               spatial_only = FALSE,
               silent = TRUE) 
#> Warning: The model may not have converged. Maximum final gradient:
#> 1.37782596110873.

# Run extra optimization steps to help convergence:
moxyb <- run_extra_optimization(moxy, nlminb_loops = 2, newton_steps = 2)
max(moxy$gradients)
#> [1] 1.377826
max(moxyb$gradients)
#> [1] 3.613655e-06

# Check the models
print(mcodb)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: ln_weight_g ~ ln_length_cm + sex_f + cpue_cod_st + year_f
#> Time column: "year"
#> SPDE: spde
#> Data: d
#> Family: student(link = 'identity')
#>              coef.est coef.se
#> (Intercept)     -4.43    0.02
#> ln_length_cm     2.96    0.00
#> sex_fM          -0.01    0.00
#> sex_fU          -0.09    0.00
#> cpue_cod_st      0.00    0.00
#> year_f1992       0.01    0.03
#> year_f1993       0.06    0.02
#> year_f1994       0.00    0.02
#> year_f1995      -0.01    0.02
#> year_f1996       0.09    0.02
#> year_f1997      -0.02    0.02
#> year_f1998      -0.06    0.02
#> year_f1999      -0.06    0.02
#> year_f2000      -0.03    0.02
#> year_f2001      -0.07    0.02
#> year_f2002      -0.08    0.02
#> year_f2003      -0.05    0.02
#> year_f2004      -0.11    0.02
#> year_f2005      -0.12    0.02
#> year_f2006      -0.10    0.02
#> year_f2007      -0.10    0.02
#> year_f2008      -0.13    0.02
#> year_f2009      -0.13    0.02
#> year_f2010      -0.12    0.02
#> year_f2011      -0.13    0.02
#> year_f2012      -0.13    0.02
#> year_f2013      -0.14    0.02
#> year_f2014      -0.14    0.02
#> year_f2015      -0.14    0.02
#> year_f2016      -0.14    0.02
#> year_f2017      -0.14    0.02
#> year_f2018      -0.16    0.02
#> year_f2019      -0.13    0.02
#> 
#> Matern range parameter: 0.63
#> Dispersion parameter: 0.08
#> Spatial SD (sigma_O): 0.09
#> Spatiotemporal SD (sigma_E): 0.06
#> Spatiotemporal AR1 correlation (rho): 0.16
#> ML criterion at convergence: -84853.842
print(moxyb)
#> Spatiotemporal model fit by ML ['sdmTMB']
#> Formula: ln_weight_g ~ ln_length_cm + sex_f + mean_DOXY_st + year_f
#> Time column: "year"
#> SPDE: spde
#> Data: d
#> Family: student(link = 'identity')
#>              coef.est coef.se
#> (Intercept)     -4.44    0.02
#> ln_length_cm     2.96    0.00
#> sex_fM          -0.01    0.00
#> sex_fU          -0.09    0.00
#> mean_DOXY_st     0.01    0.00
#> year_f1992       0.01    0.03
#> year_f1993       0.06    0.02
#> year_f1994       0.00    0.02
#> year_f1995      -0.01    0.02
#> year_f1996       0.09    0.02
#> year_f1997      -0.02    0.02
#> year_f1998      -0.06    0.02
#> year_f1999      -0.06    0.02
#> year_f2000      -0.03    0.02
#> year_f2001      -0.06    0.02
#> year_f2002      -0.08    0.02
#> year_f2003      -0.05    0.02
#> year_f2004      -0.11    0.02
#> year_f2005      -0.12    0.02
#> year_f2006      -0.10    0.02
#> year_f2007      -0.10    0.02
#> year_f2008      -0.13    0.02
#> year_f2009      -0.13    0.02
#> year_f2010      -0.12    0.02
#> year_f2011      -0.13    0.02
#> year_f2012      -0.13    0.02
#> year_f2013      -0.14    0.02
#> year_f2014      -0.14    0.02
#> year_f2015      -0.14    0.02
#> year_f2016      -0.13    0.02
#> year_f2017      -0.14    0.02
#> year_f2018      -0.16    0.02
#> year_f2019      -0.13    0.02
#> 
#> Matern range parameter: 0.62
#> Dispersion parameter: 0.08
#> Spatial SD (sigma_O): 0.07
#> Spatiotemporal SD (sigma_E): 0.06
#> Spatiotemporal AR1 correlation (rho): 0.16
#> ML criterion at convergence: -84881.355

Check the new fixed effect estimates and their confidence interval:

# Look at the new parameter (cod)
# SA: use mcod$sdreport !
# ML: Will do, thanks!
cod_est <- as.list(mcodb$sd_report, "Estimate")
cod_se <- as.list(mcodb$sd_report, "Std. Error")
cod_est$b_j[5] # 5 is index for cod parameter
#> [1] -0.0003162916
cod_est$b_j[5] + c(-2, 2) * cod_se$b_j[5] # C.I
#> [1] -0.001718512  0.001085928

oxy_est <- as.list(moxyb$sd_report, "Estimate")
oxy_se <- as.list(moxyb$sd_report, "Std. Error")
oxy_est$b_j[5] # 6 is index for oxygen
#> [1] 0.007482809
oxy_est$b_j[5] + c(-2, 2) * oxy_se$b_j[5]
#> [1] 0.005475342 0.009490275

Now compare the models using AIC:

# Compare the models with AIC (unsure actually if this is correct for a sdmTMB model...)
# SA: as long as only the fixed effects change and reml = FALSE, then yes
# SA: if random effects change and reml = TRUE, then yes
# SA: all conditional on whether AIC is a good thing to use for this type of model... but people do
# ML: I don't use restricted maximum likelihood and only fixed effects change, so I'll continue with this
aic_m1b <- extractAIC(m1b)
aic_mcodb <- extractAIC(mcodb)
aic_moxyb <- extractAIC(moxyb)

aic_m1b
#> [1]      37.0 -169633.5
aic_mcodb
#> [1]      38.0 -169631.7
aic_moxyb
#> [1]      38.0 -169686.7

The model with oxygen has the smallest AIC, and the coefficient for oxygen seems significant (unlike the cod coefficient).

Now let’s look more closely at the our estimates. If I understand Thorson (2015) correctly:

“This implies that \(\gamma X\) (the covariate times its coefficient) has a standard deviation of \(\gamma\) such that coefficients can be interpreted via comparison with the standard deviation of spatial, temporal and spatiotemporal variation, as well as that of residual variation.”

I can now compare the coefficients of the additional covariates (since they are standardized to have mean 0 and variance of 1) with the standard deviation of the spatial, spatial trend and spatiotemporal effects, i.e. \(\sigma_E\) and \(\sigma_A\) in Thorson (2015) for the spatial and spatiotemporal effects, respectively (Eqns. 6b-7). These terms are the square roots of the marginal variances of the random fields, i.e. \(\sigma_E^2\) and \(\sigma_A^2\) (denoted \(\tau^2\) above).

In sdmTMB I think \(\sigma_E\) and \(\sigma_A\) above correspond to Spatiotemporal SD (sigma_E) and Spatial SD (sigma_O) seen in print(model). I think the non-rounded values can be extracted from: moxyb$sd_report$value # SA: correct, or as a named list (note the last argument):

Ok, so for the sake of model comparison, I will now extract the estimates of 1) spatial standard deviation 2) from all models, and the coefficients of the additional covariates, together with their standard errors. The point is to compare the relative magnitude of the covariates and spatial and spatiotemporal variation, as well as comparing how the latter change with the inclusion of covariates.

# Extract coefficients from cod model
cod_est <- as.list(mcodb$sd_report, "Estimate", report = TRUE)
cod_se <- as.list(mcodb$sd_report, "Std. Error", report = TRUE)
b_cod_est <- as.list(mcodb$sd_report, "Estimate")
b_cod_se <- as.list(mcodb$sd_report, "Std. Error")

# Extract coefficients from oxygen model
oxy_est <- as.list(moxyb$sd_report, "Estimate", report = TRUE)
oxy_se <- as.list(moxyb$sd_report, "Std. Error", report = TRUE)
b_oxy_est <- as.list(moxyb$sd_report, "Estimate")
b_oxy_se <- as.list(moxyb$sd_report, "Std. Error")

# Extract coefficients from default model
def_est <- as.list(m1b$sd_report, "Estimate", report = TRUE)
def_se <- as.list(m1b$sd_report, "Std. Error", report = TRUE)

# Here are the estimates from the AIC-selected model (i.e. with cod) put in a data frame
cod_df <- data.frame(coef = c("sigma_O", "sigma_E", "cod", "oxygen"),
                     est = c(cod_est$sigma_O, cod_est$sigma_E, b_cod_est$b_j[5], NA), 
                     se = c(cod_se$sigma_O, cod_se$sigma_E, b_cod_se$b_j[5], NA),
                     model = "cod model")

# Here are the estimates from the oxygen model put in a data frame
oxy_df <- data.frame(coef = c("sigma_O", "sigma_E", "cod", "oxygen"),
                     est = c(oxy_est$sigma_O, oxy_est$sigma_E, NA, b_oxy_est$b_j[5]), 
                     se = c(oxy_se$sigma_O, oxy_se$sigma_E, NA, b_oxy_se$b_j[5]),
                     model = "oxygen (lowest AIC)")

# Here are the estimates from the default model put in a data frame
def_df <- data.frame(coef = c("sigma_O", "sigma_E", "cod", "oxygen"),
                     est = c(def_est$sigma_O, def_est$sigma_E, NA, NA), 
                     se = c(def_se$sigma_O, def_se$sigma_E, NA, NA),
                     model = "default model")

coef_df <- rbind(cod_df, oxy_df, def_df)
  
# Finally calculate C.I's
coef_df$lwr <- coef_df$est - 1.96*coef_df$se
coef_df$upr <- coef_df$est + 1.96*coef_df$se

dodge <- position_dodge(width = 0.3)

ggplot(coef_df, aes(coef, est, color = model, group = model)) + 
  geom_point(size = 2, position = dodge) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), position = dodge, width = 0.2) + 
  scale_color_brewer(palette = "Dark2")
#> Warning: Removed 4 rows containing missing values (geom_point).

I interpret this as follows: the spatial, spatial trend and spatiotemporal variation is not affected by additional covariates, apart from oxygen explaining some spatial variation that is constant in time. This makes sense since oxygen concentration largely follows depth and does not show strong spatiotemporal variation. I also find the magnitude of spatial and spatiotemporal variation to be considerably larger than that of the additional covariates.

For the sake of comparison, I can also produce a map to look at the differences there. Here I’m using the cod model.

# Add in a fixed covariate here
pred_grid_cod <- pred_grid
pred_grid_cod$mean_DOXY_st <- 0 # Mean since standardized

pcod <- predict(moxyb, newdata = pred_grid_cod)

# Replace too-deep predictions with NA
pcod <- pcod %>% mutate(est2 = ifelse(depth < -130, NA, est))

ggplot(pcod, aes(X, Y, fill = est2)) +
  geom_raster() +
  facet_wrap(~year, ncol = 6) +
  scale_fill_viridis(option = "magma", 
                     name = "log(condition factor)") + 
  geom_sf(data = world, inherit.aes = F, size = 0.2) +
  coord_sf(xlim = c(12.5, 21.5), ylim = c(54, 58)) +
  ggtitle("Prediction (random + fixed) with covariates")

To do

SA: it is a random field that represents a slope over the time variable; yes it is confusingly internally still called omega_s_slope or something similar. I should fix that. If TRUE, omega_s can be thought of as an intercept and omega_s_slope (zeta_s) the slope. We have a paper in review on it I could share if you want. Sorry it’s badly documented… although I think there is a vignette on it.

SA: extracted correctly, yes, although see my suggested list version

SA: as for that Thorson 2015 paragraph, hmm, perhaps yes, although I hadn’t thought of it that way before. I’ve more commonly seen assessment of covariates via the magnitude (ecological relevance), confidence interval, and/or AIC, but yes, that could be a useful scale comparison. That assumes you have scaled your predictor to have variance = 1, which might make less sense if it’s a log variable.

Updated

References

Anderson, S.C., Keppel, E.A., Edwards, A.M. 2019. A reproducible data synopsis for over 100 species of British Columbia groundfsh. DFO Can. Sci. Advis. Sec. Res. Doc. 2019/041. vii + 321 p.

Casini, M., Käll, F., Hansson, M., Plikshs, M., Baranova, T., Karlsson, O., Lundström, K., Neuenfeldt, S., Gårdmark, A. and Hjelm, J., 2016. Hypoxic areas, density-dependence and food limitation drive the body condition of a heavily exploited marine fish predator. Royal Society open science, 3(10), p.160416.

Froese, R., Thorson, J.T. and Reyes Jr, R.B., 2014. A Bayesian approach for estimating length‐weight relationships in fishes. Journal of Applied Ichthyology, 30(1), pp.78-85.

Grüss, A., Gao, J., Thorson, J.T., Rooper, C.N., Thompson, G., Boldt, J.L. and Lauth, R., 2020. Estimating synchronous changes in condition and density in eastern Bering Sea fishes. Marine Ecology Progress Series, 635, pp.169-185.

Gårdmark, A., Casini, M., Huss, M., van Leeuwen, A., Hjelm, J., Persson, L. and de Roos, A.M., 2015. Regime shifts in exploited marine food webs: detecting mechanisms underlying alternative stable states using size-structured community dynamics theory. Philosophical Transactions of the Royal Society B: Biological Sciences, 370(1659), p.20130262.

Neuenfeldt, S., Bartolino, V., Orio, A., Andersen, K.H., Andersen, N.G., Niiranen, S., Bergström, U., Ustups, D., Kulatska, N. and Casini, M., 2020. Feeding and growth of Atlantic cod (Gadus morhua L.) in the eastern Baltic Sea under environmental change. ICES Journal of Marine Science, 77(2), pp.624-632.

Orio, A., Bergström, U., Florin, A.B., Lehmann, A., Šics, I. and Casini, M., 2019. Spatial contraction of demersal fish populations in a large marine ecosystem. Journal of Biogeography, 46(3), pp.633-645.

Orio, A., Bergström, U., Florin, A-B., Sics, I. and Casini, M., 2020. Long-term changes in spatial overlap between interacting cod and flounder in the Baltic Sea. Hydrobiologia, 847(11), pp.2541-2553.

Svedäng, H. and Hornborg, S., 2014. Selective fishing induces density-dependent growth. Nature communications, 5(1), pp.1-6.

Thorson, J.T., 2015. Spatio-temporal variation in fish condition is not consistently explained by density, temperature, or season for California Current groundfishes. Marine Ecology Progress Series, 526, pp.101-112.